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We report on the behavior of the eigenvalue distribution of the Dirac operator in (2+l)-flavor QCD 
at finite temperature, using the HISQ action. We calculate the eigenvalue density at several values 
of the temperature close to the pseudocritical temperature. For this study we use gauge field 
configurations generated on lattices of size 32^ x 8 with two light quark masses corresponding to 
pion masses of about 160 and 115 MeV. We find that the eigenvalue density below Tc receives 
large contributions from near-zero modes which become smaller as the temperature increases or 
the light quark mass decreases. Moreover we find no clear evidence for a gap in the eigenvalue 
density up to lATc. We also analyze the eigenvalue density near Tc where it appears to show a 
power-law behavior consistent with what is expected in the critical region near the second order 
chiral symmetry restoring phase transition in the massless limit. 
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1. Introduction 

The chiral phase transition is one of the most important features of quantum chromodynamics 
(QCD) in our understanding of the properties of strongly interacting matter in a hot medium. QCD 
with A//-flavors of massless quarks has an SUiiNf) x SUR(Nf) chiral symmetry, which is sponta- 
neously broken in the vacuum, i.e., at zero temperature. The vanishing of the chiral condensate 
(y/y/), an order parameter of the chiral phase transition, signals the restoration of chiral symmetry 
at high temperatures. 

It was shown by Banks and Casher ^ that the eigenvalue density of the Dirac operator p(A) 
is related to (xf^yf) by \{\l^y/) \ = ^p(O) in the limits of infinite volume and vanishing quark mass. 
Therefore low-lying eigenvalues of the Dirac operator play a key role in the breaking and restoration 
of chiral symmetry. In particular, p (0) is nonzero in the chirally broken phase and p (0) vanishes 
when the chiral symmetry is restored. Whether p(A) develops a gap around A = above the 
critical temperature Tc is another important subject since it is related to restoration of the C/a(1) 
symmetry. The C/a(1) symmetry is explicitly broken by the anomaly which exists independently 
of temperature. Thus the C/a(1) symmetry is expected to remain broken above Tc. If the difference 
of isovector susceptibilities co = Xp - Xs [^l, where Xi = / d^^{jK^)JH^)) ^ ^i^)^ jp{^) = 
x[r{x)T^iY5\l/{x) and = v/(x)t^i//(x), is considered, a relation \ 



where m is a quark mass, can be obtained in a manner similar to the Banks-Casher relation. There- 
fore, if p (A) has a gap around A = 0, then = in the chiral limit, which means Xp = Xs and that 
the C/a(1) symmetry is effectively restored. 

If the chiral phase transition occurs continuously and p(A) changes smoothly as the tempera- 
ture increases, a power-law behavior of p(A) such as p(A) ^ A^ becomes important because a is 
related to critical exponents of the phase transition. It has been suggested that for Nf = 2, there is 
a second order phase transition [||] belonging to the same universality class as 3-dimensional 0(4) 
spin models. In fact, for staggered quarks, although only an 0(2) symmetry is preserved rather 
than 0(4), at least 0{N) scaling has been suggested in an analysis of the magnetic equation of 
state 0, ^. Similarly, studies performed with Wilson quarks indicated 0(4) scaling [^0]. Thus an 
interesting question is whether the power-law of p(A) at the critical temperature is consistent with 
what is expected from the scaling behavior for the 0(2) or 0(4) universality class. 

In this study, we calculate low-lying eigenvalues of the Dirac operator at several values of 
the temperature close to the pseudocritical temperature by using (2+l)-flavor dynamical quarks 
with the highly improved staggered quark (HISQ) action. In the following sections, we show the 
temperature and light quark mass dependence of p(A) and its critical behavior. 

2. Simulation setup 

Our simulations were performed on 32^ x 8 lattices with the tree level improved gauge action 
and the HISQ action which reduces the effects of taste symmetry violations and accordingly the 

^Here and in the following we ignore any explicit dependence of p(A) on the light sea quark mass m, i.e. we write 
p(A) = p(A,m). This will eventually become important also for a discussion of the critical behavior. 
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Table 1: The parameters of the numerical simulations, i.e. the gauge coupling j8, the lattice spacing a, tem- 
perature T and the strange quark mass are summarized. The number of configurations for measurements 
is also listed. 

cutoff dependence better than the other staggered fermion formulations in use [|[]. Part of our 
gauge configurations were generated by the HotQCD collaboration The lattice spacing was 
determined by measuring the static quark anti-quark potential. The strange quark mass was set 
to its physical value and the light quark mass m/ was fixed to m^/20 and m^/40 corresponding to 
lightest (Goldstone) pions of about 160 and 115 MeV, respectively. The pseudocritical temperatures 
for mi/rus = 1/20 and 1/40 has been estimated at 162.9(1.8) and 157(3) MeV, respectively, from 
the peak in the chiral susceptibility. The details of the determination of the lattice spacing, the 
strange quark mass and the pseudocritical temperature has been discussed in Ref. [^. Each 10th 
trajectory was chosen for measurements after skipping at least 500 trajectories for thermalization. 
Statistical errors were estimated by the jackknife method. Our simulation parameters and statistics 
are summarized in Table [T]. 

The staggered Dirac operator^ is anti-hermitian and it has purely imaginary eigenvalues iX. 
Because of the remnant chiral symmetry, the eigenvalues always appear in complex conjugate pairs. 
We calculated the lowest 100 positive eigenvalues Ay^ (^ = 1, 2, • • • , 100) of il/>, and then evaluated 
the eigenvalue density defined by 

p(A)^l(^5(A-A,)), (2.1) 

^ k 

where V is the four- volume. Here p(A) is normalized so that / JAp(A) = # eigenvalues/V . Nu- 
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Figure 1: Temperature and light quark mass dependence of the eigenvalue density. Plus and cross symbols 
indicate p(A) for m//m^ = 1/20 and 1/40, respectively. The value of the m/'s are shown by vertical dashed 
and dotted lines. 



merically, p(A) was computed by binning of eigenvalues in small intervals for each configuration. 
In this study, we chose 0.0005 as the size of a bin for all of the temperatures and light quark masses. 

3. Temperature and light quark mass dependence of p(A) 

Figure [l] shows the temperature and light quark mass dependence of p(A). We found a small 
quark mass dependence for eigenvalues A small compared to m/, indicated by vertical dashed and 
dotted lines in the plots. The bulk behavior of p (A), on the other hand, is not sensitive to the quark 
mass. Moreover, p(0) becomes smaller as the temperature increases and/or the light quark mass 
decreases, and eventually goes to zero, consistent with what is expected from the restoration of 
chiral symmetry, i.e., the vanishing of the chiral order parameter above T^. 

The eigenvalue distribution p(A) at temperatures clearly above Tc is shown in Fig. ^ The 
distribution of eigenvalues above T = 239.7 MeV is qualitatively different from that below that 
temperature. Namely, p(A) has only a tail that approaches the origin. As mentioned in Sec. |l|, the 
important point here is whether p(A) has a gap around A = 0. From the r.h.s of Fig. a lack of 
eigenvalues around the origin can be seen above T = 188.7 MeV. However, a more careful analysis, 
including an investigation of the (spatial) volume dependence, is needed to distinguish a real gap 
from lattice artifacts and finite volume effects. 

4. Critical behavior of p(A) 

In the vicinity of a critical point, the behavior of an order parameter M is controlled by a 
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Figure 2: The eigenvalue density above Tc. Plus, cross, asterisk, square, circle, triangle, and downward 
triangle symbols indicate p(A) at T =173.0, 177.7, 188.7, 210.6, 239.7, 275.9, and 331.6 MeV, respectively. 
The right figure shows a logarithmic plot of the same data for the small A region. All of the data points for 
r = 331.6 MeV are outside of the region shown. 



universal scaling function /g(z) as 

M = /z1/Vg(z), (4.1) 

with z = t/h^l^^, where h and t are scaling variables corresponding to a symmetry breaking field 
and temperature, respectively, and j8 and 5 are critical exponents. In QCD (y/i//) and the (light) 
quark mass m are regarded as M and /z, respectively. Thus one has the relation 

{WW)-m''^fG{z) . (4.2) 

On the other hand, in the infinite volume limit, {'^f^) can be obtained from the eigenvalues of 
the Dirac operators as 

2mp{X) 



(V^V^> = -MA^. (4.3) 







Assuming p (A) ^ AA^, Eq. (|4.3| ) can be rewritten in the limit m ^ as 

2AA« 



(V/v^) = -m«/ dlj:^ (4.4) 
Jo A^ + 1 



with A = A/m. Thus, by comparing to a = 1/5 would be expected at Tc in the chiral 
limit and a should have a value close to 1/5 for a small enough quark mass and near Tc. 

To test this expectation, we fit the eigenvalue density around Tc to the Ansatz p(A) = AA^. 
Here we set the fit range as [0, Amax]- Since the part of p(A) with large A is suppressed due to us 
having calculated only a fixed number of low-lying eigenvalues per configuration, the largest A in 
the region without such a suppression effect is chosen as Amax- 

Figure ^ shows the temperature dependence of the fit parameters a and A. a increases mono- 
tonically as the temperature increases and it has a value close to 1/5 for either the 0(2) or 0(4) 
universality class^ at a temperature not more than 10 MeV below the pseudocritical temperature 
for both mi/rus = 1/20 and 1/40. Since we expect that a = 1/5 at only in the chiral limit, the 



' 1/5 for the 0(2) and 0(4) universality classes are too similar to be distinguishable within our numerical accuracy. 
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Figure 3: Temperature dependence of the fit parameters a (left) and A (right). Plus and cross symbols 
correspond to mi/rus = 1/20 and 1/40 results, respectively. In the left figure, the inverse of the critical 
exponent 5 for the 0(2) [ p^ ] and 0(4) [|ll| universality classes is indicated by the horizontal dashed and 
dotted lines, respectively. Note that those two horizontal lines are too close to each other to be distinguished. 
The pseudocritical temperature for each light quark mass is shown by a colored band. 




T [MeV] T [MeV] 

Figure 4: ;^^/dof for our fit analysis with the Ansatz p(A) = AX^ at each temperature (left). ;^^/dof in the 
case that a is fixed to 1/5 for the 0(2) universality class is shown in the right figure — note the different 
y-axis scale. The pseudocritical temperature for each light quark mass is shown by a colored band. 



fact that this occurs at a somewhat smaller temperature should be due to the finite m/ used. The 
deviation becomes smaller as m/ decreases. We also note that for finite m/, namely in the crossover 
region, the (pseudo-)critical temperature depends on the quantity which is used to determine it. A 
shows a monotonically increasing behavior, too, but it looks insensitive to temperature below T^. 

We show the ;^^/dof of our fit analysis in Fig. ^. We find a minimum of ;f^/dof with a value 
close to 1 for both of the quark masses near their pseudocritical temperatures. Especially, if we fix 
a to 1/5 for the 0(2) universality class, /Aof is quite large when the temperature is far from 
the pseudocritical one^. This means that our fit analysis and Ansatz work well only around the 
pseudocritical temperature, which is consistent with what is expected from the scaling behavior of 
the order parameter discussed above. 

^The plot using the 0(4) value of 1/5 would be indistinguishable. 
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5. Conclusions 

We studied the behavior of eigenvalue distributions below and above Tc with (2+l)-flavor dy- 
namical quarks using the HISQ action. We find that the eigenvalue density around A = decreases 
as the temperature increases and/or the light quark mass decreases. This is consistent with what is 
expected from the behavior of the chiral order parameter. Moreover, we find that the eigenvalue 
density above Tc has a tail that approaches the origin and that there is no clear evidence for a gap 
in the eigenvalue density up to 177.7 MeV, about I.IT^. 

We also investigated the critical behavior around Tc by fitting the eigenvalue density to an 
Ansatz p(A) = AA^, which is the form expected from the scaling behavior of the order parameter 
around Tc. We find that the fit analysis works well around the pseudocritical temperature and that 
the value of a there is compatible with 1/5, where 5 is the critical exponent for the 0{2) or 0(4) 
universality class which governs the behavior of the order parameter ({xj/Xj/)) at Tc. 
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